Heart rate variability changes with respect to time and exercise intensity during heart-rate-controlled steady-state treadmill running

The aim of this work was to investigate the time and exercise intensity dependence of heart rate variability (HRV). Time-dependent, cardiovascular-drift-related increases in heart rate (HR) were inhibited by enforcing a constant heart rate throughout the exercise with a feedback control system. Thirty-two healthy adults performed HR-stabilised treadmill running exercise at two distinct exercise intensity levels. Standard time and frequency domain HRV metrics were computed and served as outcomes. Significant decreases were detected in 8 of the 14 outcomes for the time dependence analysis and in 6 of the 7 outcomes for the exercise intensity dependence analysis (excluding the experimental speed-signal frequency analysis). Furthermore, metrics that have been reported to reach an intensity-dependent near-zero minimum rapidly (usually at moderate intensity) were found to be near constant over time and only barely decreased with intensity. Taken together, these results highlight that HRV generally decreases with time and with exercise intensity. The intensity-related reductions were found to be greater in value and significance compared to the time-related reductions. Additionally, the results indicate that decreases in HRV metrics with time or exercise intensity are only detectable as long as their metric-specific near-zero minimum has not yet been reached.

www.nature.com/scientificreports/ HRV frequency bands including very-low (VLF) and ultra-low frequency (ULF) 11 . Consistent with previous reports 10 , HRV was found to decrease with increasing exercise intensity, and HRV was observed to decrease over time. But the confounding effects of cardiovascular drift on the time-dependent HRV decrease were still at play, and the study concluded that it remains to be clarified whether these changes are due to time or due to increases in HR related to cardiovascular drift 11 . In the work proposed here, to overcome the drift-related limitations, a feedback control system was employed to keep the heart rate constant throughout the duration of treadmill exercise 12,13 . The technical feasibility of using HR control for the investigation of time-related changes in HRV was proven previously in a pilot study with eight participants 14 ; while the outcome measures showed an overall tendency to decrease over time, the decrease was significant ( p < 0.05 ) for only seven of the 10 HRV metrics, pointing to the study being statistically underpowered. Further limitations identified in the pilot study, and which are eliminated in the present work, included the employment of a feedback design with low-pass characteristics, thus attenuating information content in the LF and HF bands, and the reconstruction of RR intervals from the HR signal, rather than from direct recordings of raw RR intervals.
By forcing the heart rate to stay constant for the duration of the treadmill exercise, confounding HR-related influences can be removed, and the unobstructed HRV analysis becomes viable. A key challenge with this new HR control approach for HRV analysis is the effect the control loop's compensating actions can have on specific frequency ranges, potentially affecting frequency-related HRV metrics. Following a suggestion made in our pilot work 14 , a distinct feedback control structure was implemented here to address this limitation by using an inputsensitivity-shaping approach to obtain a uniformly flat frequency response across all frequencies of the treadmill speed signal. Therefore, the treadmill speed signal might be able to act, instead of or as well as the recorded RR intervals, as an indirect source for the frequency-domain HRV analysis, where frequency-specific control loop compensation effects are eliminated. The viability of using speed in this way is explored in the sequel.
The fundamental idea of implementing feedback control of HR to stabilise HR for time-dependency analysis of HRV, as well as the concept of using the treadmill speed signal as a proxy for the RR interval signal to perform an indirect frequency-domain HRV analysis, is new. Aside from our pilot study 14 , this has not been investigated in any previous research.
The aim of this work was to investigate the time and exercise intensity dependence of HRV during steady-state treadmill running while using feedback control to prevent HR drift. Based on the previous findings of Michael et al. 10 and Hunt et al. 11 , we hypothesised that HRV can be expected to decrease with increasing exercise intensity and throughout the duration of the exercise.

Methods
Participants. Thirty-two healthy, regularly exercising (three times per week and at least 30 min per session) adults were included in the study (Table 1; details of an a priori sample size estimate are given later, "Sample size estimate"). Recruitment was carried out by convenience sampling within the Department of Engineering and Information Technology of Bern University of Applied Sciences in Burgdorf and the University of Bern. Of the 32 participants, 29 were male and 3 were female. Smokers and persons with prior history of cardiovascular or respiratory disease or current musculoskeletal complaints or injuries were excluded. Before each test, participants were required to avoid strenuous exercise (24 h), caffeine (12 h), and heavy meals (4 h). The study was approved by the local ethics committee (Ethics Committee of the Swiss Canton of Bern, Ref. 2021-00889), and the participants provided written informed consent prior to participation.
Feedback control design. The control structure was set up as a generic negative feedback control system ( Fig. 1) with feedback compensator C(s) and nominal plant P o (s). It is known that the plant parameters k and τ vary widely between different people, yet it has been clearly demonstrated that a constant-coefficient, linear, time-invariant compensator based on an average nominal model can deliver stable and accurate HR control performance when applied to a wide range of participants [12][13][14] , that is to say, robust control is obtained as a consequence of the fundamental ability of feedback to reduce plant uncertainty 15 .
The nominal plant parameters employed here were set to average values obtained from a total of 73 participants from two separate model identification studies, namely k = 24.88 and τ = 59.28 . These values were derived by combining parameters obtained empirically in two previous studies, more specifically k = 26.2, τ = 62.5 with n = 25 16 and k = 24.2, τ = 57.6 with n = 48 17 ; the respective sample sizes were taken into account.
The linear feedback compensator, generally described in rational form as C(s) = G(s)/H(s) , was intentionally constructed as a merely proper first-order transfer function with an integral term, viz.
where n g = 1 and n h = 1 are the degrees of polynomials G and H. As detailed below, this combination of nominal plant Eq. (1) and compensator Eq. (2) results in a closed-loop characteristic polynomial of degree two and, by virtue of the compensator's two free parameter g 1 and g 0 , allows arbitrary placement of the two closed-loop poles in the complex plane.
The choice of a proper transfer-function for C is to ensure that the input sensitivity function U o (Eq. (5), below) is also proper and thus remains finite over all frequencies, in line with our design goal to make |U o | constant for all frequencies (a strictly proper C would make U o strictly proper, whence |U o | would roll off towards zero as frequency increases).
The closed-loop sensitivity function S o , complementary sensitivity function T o and the input sensitivity function U o can be deduced from the control structure 15 . S o describes the transfer function from disturbance d to the controlled output y ( d → y ), T o from the reference signal r to the controlled output y ( r → y ), and U o from the reference signal r and the disturbance d to the treadmill speed signal u ( d → u and r → u): The closed-loop characteristic polynomial can be identified as The feedback design goal is to shape the input-sensitivity magnitude to be constant across all frequencies. The rationale for this goal is that the treadmill speed command signal u, which is linked to HRV disturbance term d through U o , might potentially be used for frequency-domain HRV analysis across the whole frequency range.
, and we proceed by taking a cancellation approach to simplify this down to a constant value. As a first step, all plant poles are cancelled by constraining the compensator numerator polynomial G to include the known factor A o and a remaining, unknown factor G ′ by writing G = A o G ′ , which results in Since the plant is assumed strictly proper, n b < n a , and the compensator proper, n g = n h , it follows that the order of the characteristic polynomial, Eq. (6), is n φ = n a + n h . Since a unique solution of Eq. (6) for G and H requires that n φ be equal to the number of free parameters in G and H, namely n h − 1 + n g + 1 = n h + n g , it follows that n a + n h = n h + n g ⇒ n g = n a (giving also n h = n a , since n g = n h ). This in turn implies that, since by construction G = A o G ′ , the polynomial G ′ must be of degree zero, i.e. G ′ is a constant, which we denote g ′ 0 . This allows U o to be written, using Eq. (7), as The compensator parameters in G(s) = g 1 s + g 0 can be identified by noting G = g ′ 0 A o = 1 k (s + 1 τ ) to give g 1 = 1/k and g 0 = 1/(kτ ) , wherefore and it is seen that the compensator parameters depend only on the plant parameters. Using the nominal values k = 24.88 and τ = 59.28 , the specific compensator used in this study is calculated as The sensitivity functions S o , T o and U o for this compensator with the nominal plant parameters can be computed using Eqs. (12)- (14) and are plotted in Fig. 2. In particular, it can be seen that |U o (jω)| (blue line) maintains a constant value of 1/k over all frequencies ( k = 24.88 and 20 log 10 (1/24.88) = −28 dB) and that the bandwidth of both S o and T o is the same as the open-loop plant bandwidth ( τ = 59.28 and 1/(59.28 × 2 × π) = 0.0027 Hz). We remark, for the purposes of later discussion, that this bandwidth is just below the upper bound of the ULF frequency band which lies at 0.0033 Hz.
Experimental design and testing protocol. Each participant performed two treadmill running tests at different exercise intensities (lower intensity level 1 [ EIL 1 ] and higher intensity level 2 [ EIL 2 ]), and each with a duration of 35 min. The intensities for the two tests were set for each participant individually. Each test was carried out on a separate day with at least 48 h between tests. The recording of the two tests ( EIL 1 and EIL 2 ) was counterbalanced by randomizing the order of employed intensity in order to avoid any order-of-presentation effects. A 10-min warm up was used to assess the HR level and speed at a given baseline exertion level. Participants were asked to choose a running speed equivalent to an exertion level of 13 (i.e. somewhat hard) on the Borg rating of perceived exertion (RPE) scale 18 . The HR recorded during the warm up ( HR Warmup ) was then averaged and used to calculate the two exercise intensities EIL 1 and EIL 2 . The reference HR for each EIL was calculated as a symmetric deviation from HR Warmup by 3% of the maximal age-related heart rate ( HR max = 220 − age ) as Outcome measures. Analysis. Time dependence analysis was performed by calculating and comparing HRV metrics for two consecutive evaluation windows ( w 1 and w 2 ) with equal duration of 15 min: window w 1 was defined from t = 5 min to t = 20 min and window w 2 from t = 20 min to t = 35 min . For the exercise intensity dependence analysis, HRV metrics were computed for the total duration (denoted w 1 w 2 ) of EIL 1 and then compared to the HRV metrics computed for the total duration ( w 1 w 2 ) of EIL 2 . RR interval outliers or ectopic heartbeats 19 were identified and removed using an impulse rejection filter 20 with a threshold value of 7 and a window length calculated to span about 25 s.

Time domain HRV metrics.
In the time domain, the root-mean-square of successive differences (RMSSD) and the standard deviation of the normal-to-normal (NN) intervals (SDNN) were evaluated and served as primary endpoints. SDNN and RMSSD are defined as follows: where NN i is the ith recorded interval and NN is the mean value of all NN intervals. N is the total number of NN intervals.
Frequency domain HRV metrics. In the frequency domain, the average power contained in the frequency bands ULF, VLF, LF, HF, and the total frequency range (TP, Total Power), was computed for the RR signal and the treadmill speed signal using the Lomb-Scargle method for spectral density estimation.
Statistical analysis. Sample size estimate. An a priori statistical power calculation using observed effect size and dispersion from a pilot study with a sample of n = 8 participants 14,21 was performed to obtain a sample-size estimate. Calculations used a significance level of α = 0.05 and a required statistical power of 1 − β = 0.9 . The statistical power calculation was based on the total power of the HRV spectrum (TP, total power spectral density estimate). The pilot study found a difference in paired samples for the log-transformed data of 0.133 ± 0.236 (mean ± standard deviation). The power calculation revealed a required sample size of n = 29 . To cater for the possibility of participants dropping out, a ∼10% contingency buffer (+3 participants) was added, resulting in a final sample size of n = 32.
Analysis. Based on our hypothesis that HRV can be expected to decrease with increasing exercise intensity and throughout the duration of the exercise, one-sided t-tests with significance level α = 0.05 were performed for www.nature.com/scientificreports/ all HRV metrics (time and frequency domain) to assess time and exercise intensity dependence. For the time dependence analysis, HRV metrics computed for w 1 were compared to the metrics for w 2 (for both intensity levels EIL 1 and EIL 2 ). For the exercise intensity dependence analysis, HRV metrics for the total duration of EIL 1 were compared to those for the total duration of EIL 2 . All metrics were log-transformed to reduce skew and make the data conform more closely to the normal distribution. A Shapiro-Wilk test was performed to assess normality, and no data was found to be significantly different from normal.

Equipment and data collection. All tests were performed on a computer-controlled treadmill (model
Pulsar, h/p/cosmos Sports & Medical GmbH, Germany). Treadmill speed was set by a computer running the HR feedback controller within the Simulink Desktop Real-Time environment (The MathWorks, Inc., USA) and communicating with the treadmill over the coscom v3 interface protocol. Heart rate and raw RR intervals were recorded with a chest-strap-mounted sensor (H10, Polar Electro Oy, Finland) transmitted to an ESP32 development board (Espressif Systems, China) over Bluetooth low energy and sent via serial communication to the control application running on the PC. A Polar V800 wristwatch was employed as a backup method for saving the heart rate and RR intervals. The control application was set up to work with heart rate values transmitted in units of beats per minute (bpm). The V800 and the ESP32 development board saved the more accurate RR intervals, recorded with millisecond resolution, for later analysis.

Results
Sixty-two of the 64 measurements were recorded successfully (one participant dropped out of the study due to muscle pain). Fifty-five of the 62 recorded data sets were able to be used for analysis: three measurements had to be excluded due to insufficient HR control performance, another three due to HR abnormalities and one due to faulty sensor data.
For illustration, a sample data set for a single participant is provided (Fig. 3), while a complete set of categorical scatter plots showing all hypothesis testing results is provided for time domain (Fig. 4) and frequency domain (Fig. 5) outcomes. Table 2 lists the corresponding numerical values.
Time domain outcomes. SDNN was found to significantly decrease over time at EIL 1 ( p = 0.007 ) and EIL 2 ( p < 0.001 ), and was observed to be lower at a higher intensity ( p < 0.001 ). RMSSD, on the other hand, revealed no significant differences, neither for the intensity comparison ( p = 0.209 ) nor the time dependence analysis at EIL 1 ( p = 0.527 ), nor the time dependence analysis at EIL 2 ( p = 0.587).
Frequency domain outcomes. The frequency domain outcomes are separated into results derived from the RR interval analysis and results derived from the speed signal analysis (Fig. 5).
For the speed signal frequency domain analysis, ULF, VLF, and TP were found to decrease over time at both EIL 1 and EIL 2 ( p, ULF Speed < 0.001 ; VLF Speed = 0.003 ; TP Speed < 0.001 for EIL 1 and p, ULF Speed = 0.009 ; VLF Speed < 0.001 ; TP Speed < 0.001 for EIL 2 ). A significant difference was found in LF at EIL 2 ( p, LF Speed = 0.002 ), and moderate evidence for a decrease was observed in LF at EIL 1 ( p, LF Speed = 0.082 ). No significant differences in HF at EIL 1 ( p, HF Speed = 0.147 ) and moderate evidence at EIL 2 ( p, HF Speed = 0.071 ) were identified. The exercise intensity dependence analysis revealed significant differences in LF and no evidence in ULF, VLF, HF, and TP ( p, ULF Speed = 0.997 ; VLF Speed = 0.434 ; LF Speed = 0.004 ; HF Speed = 0.780 ; TP Speed = 0.982).

Discussion
This work aimed to investigate the time and exercise intensity dependence of HRV during steady-state treadmill running while using feedback control to prevent HR drift: we hypothesised that HRV can be expected to decrease with increasing exercise intensity and throughout the duration of the exercise.
Analysing the time dependence of time and frequency domain HRV metrics, EIL 1 and EIL 2 revealed similar findings. SDNN, as well as ULF, VLF, LF, and TP (all except HF), computed using RR intervals and the treadmill speed signal, in the main showed strong evidence of decreasing in value over time (with moderate evidence for ULF RR at EIL 2 → p = 0.072 ; LF RR at EIL 1 → p = 0.087 and LF Speed at EIL 1 → p = 0.082 ). RMSSD and HF (RR and speed) did not significantly decrease over time. Combining these results with the tendency that most HRV measures reach an intensity-dependent near-zero minimum, and that HRV metrics believed to be associated with cardiac parasympathetic activity (i.e. RMSSD and HF) have been reported to reach that minimum more rapidly (usually at moderate intensity) 10 , gives good reason to infer that HRV decreases with time as long as the intensity-dependent near-zero minimum has not yet been reached. This conclusion is also supported by the www.nature.com/scientificreports/ Table 2. Numerical outcomes. All outcomes are log-transformed, thus dimensionless; MD mean difference (intra-group comparison: w 2 − w 1 ; inter-group comparison: EIL 2 -EIL 1 ); CI upper 95% confidence interval boundary; p value: derived from a single-sided t-test performed on the log-transformed data of the respective comparison groups. The p values were conditionally emphasized: p < 0.05 → italic else bold).     14 , where exercise performed at light intensity revealed HF and an RMSSD-related metric to decrease over time ( p = 0.053 for RMSSD proxy; p = 0.047 for HF). It would seem that with the low exercise intensity, RMSSD and HF had not yet reached a minimum and thus were able to decrease over the observation interval.
The intensity-level comparison identified a significant decrease in SDNN and in all frequency-domain HRV metrics (computed for the RR intervals) with increasing exercise intensity. RMSSD did not change, indicating the presence of an intensity-dependent near-zero minimum below EIL 1 . Despite the consistent decrease in frequency-domain HRV metrics computed for the RR intervals, the corresponding outcomes for the speed signal's frequency-domain HRV metrics differed greatly. VLF and HF stayed approximately constant, while ULF and TP increased. This finding is contrary to expectations considering the connection of the RR interval signal via the controller to the speed signal. With the control loop's flat input sensitivity function and the consistent control structure across both exercise intensities, intensity-dependent changes in outcomes were expected to be reflected in both the RR intervals and the speed signal. Based on the large increase in ULF power, we suspect this outcome to be affected by an inconsistent plant model. A deviation from the nominal plant can directly affect the actual characteristics of the sensitivity functions. The conspicuous increase in ULF leads us to believe that a decrease in the actual steady-state gain parameter k with increasing exercise intensity might have been the cause as the limit of |U o (jω)| as frequency tends to 0 is 1/k (see "Feedback control design"). A lower than nominal k would reduce the control loop's dampening effect on the impact the RR signal can have on the treadmill speed, thus leading to an overall increase in power. However, this argument may not be in agreement with a previous study 17 , where model plant parameters k and τ were reported not to be significantly dependent on exercise intensity. A second contributory factor might have been the larger speed reduction observed at higher intensities. Throughout a running exercise bout, the treadmill speed typically drifts downwards to compensate for cardiovascular-driftrelated HR changes: this downward trend in treadmill speed was observed to be greater at higher intensities.
As a consequence, we suggest that the method of using the treadmill speed signal as a proxy for the RR intervals during heart rate stabilised running exercise, in order to perform frequency-domain HRV analysis where frequency-specific control loop compensation effects are eliminated using input-sensitivity-shaping, needs further investigation to be applicable for an intensity dependence analysis. On the other hand, for the time dependence analysis, this method produced results that closely matched the trends found in frequencies primarily unaffected by the control loop (VLF, LF, and HF) of the original RR interval analysis. This suggests that the speed signal could be used as an alternative to detect trends in frequencies primarily affected by the control (namely ULF) without the influence of the control loop's compensation effects.
We acknowledge that the duration of 15 min for the analysis windows w 1 and w 2 is relatively short for the ULF power estimation. Longer analysis windows are generally more desirable but would place additional demands on the participants performing the running exercise, therefore a balance has to be found.

Conclusion
In summary, feedback control of heart rate was successfully employed to answer the question of whether timedependent HRV changes occur due to time itself or due to cardiovascular-drift-related heart rate increases. Most HRV metrics were found to decrease with time and with exercise intensity. The exercise-intensity-related reductions were generally found to be greater in value and significance compared to the time-related reductions. HRV metrics that have been reported to reach an intensity-dependent near-zero minimum rapidly (usually at moderate intensity) were found to be near constant over time and only barely decreased with intensity, indicating that decreases in HRV metrics with time or exercise intensity are only detectable as long as their metric-specific near-zero minimum has not yet been reached.

Data availability
The raw data supporting the conclusions of this article can be found in the OLOS repository (https:// doi. org/ 10. 34914/ olos: jrpnz 3eeh5 ephkl jdcai eefxqi) and will be made available by the authors on reasonable request.